clear    ;
close all;
%% ----- GRID OVER gamma -----
gamma_v = 0.1:0.1:0.9;% 0.3;% 
nG = numel(gamma_v);

% Preallocate a struct array to store everything per gamma
results(nG) = struct('gamma',[], 'ss',[], 'volat',[], ...
    'loss_peg',[], 'loss_noFXI',[], 'loss_pegsub',[],'loss_noFXI_dcp',[],'loss_suboptimalFXI',[], ...
    'lambda_peg',[], 'lambda_noFXI',[], 'lambda_pegsub',[],'lambda_noFXI_dcp',[], 'lambda_suboptimalFXI',[],...
    'CHss',[], 'Yss',[], 'CHstarss',[], 'LOss',[], ...
    'Delta_lnq',[], 'Ld_e',[],'Ld_e_dcp',[],'Delta_lnrer',[],'Delta_lnrer_dcp',[],'weight',[],'uwedge',[],'vwedge',[]);


for jj = 1:nG
    fprintf('\n=== Running gamma = %.3f (jj=%d/%d) ===\n', gamma_v(jj), jj, nG);

    %% ---------- PARAMETERS (exogenous to the fixed-point step) ----------
    Cstar = 0.5;   % to have oil export to total export ratio of 0.62 (GCC average 2012-2023, GCCSTAT) 
    Ow    = 0.506; % to get 0.70 non-oil GDP to GDP ratio (2000-2022, WB)
    Ass   = 1.0;
    gamma = gamma_v(jj);   % trade openness (looped)
    sigman = 0.10;
    alpha = 1 - 0.17;
    epsx  = 1.5;  % Camara et al. (2024) and many others
    beta  = 0.99;
    omega = 345;

    epsh  = 6;
    theta = 0.75;           % price stickiness
    Rss   = 1/beta;

    rhor  = 0.9; rhon=0.90; rhoa=0.90; rhopos=0.9; rhors=0.9; rhopf=0.9;
    eta   = (epsh-1)/epsh;
    Rstarb = 1/beta;
    gammaPi = 2.0;
    gammaYY = 0*0.05;
    Piss    = 1;

    params = [Cstar; Ow; gamma; alpha; epsx; Ass];

    %% ---------- STEADY STATE (FSOLVE) ----------
    options = optimset('Display','iter');
    [ss,Fval,exitflag] = fsolve(@(x) OEXRPolicy_ss(x, params), 8.52, options);

    % Baseline steady-state objects
    PFstarss=1; POstarss=1; PIss=1; PHss=1;

    LOss     = ss;
    CFss     = gamma * (alpha)/(1-alpha) * POstarss/(LOss)*1/PFstarss;
    PHstarss = epsx/(epsx-1) * POstarss/((1-alpha)*Ass*LOss^alpha);
    CHstarss = PHstarss^(-epsx)*Cstar;
    Oss      = PHstarss*CHstarss/POstarss - PFstarss*CFss/POstarss + Ow;
    CHss     = Ass*LOss^alpha*Oss - CHstarss;
    ess      = gamma/(1-gamma) *CHss/CFss; %ess = (1-alpha)*PIss*Ass*LOss^alpha/POstarss; 
    Yss      = Ass * LOss^alpha * Oss;
    Lss      = LOss*Oss;
    rerss    = ess;


    
    NX_nonoil_output  = ess*( PHstarss*CHstarss - PFstarss*CFss )/ Yss;
    EX_oil_output     = ess*POstarss*(Ow - Oss)                  / Yss;
    EX_nonoil_output  = ess*PHstarss*CHstarss                    / Yss;
    IM_nonoil_output  = ess*PFstarss*CFss                        / Yss;
 
    GDP           = Yss + ess*POstarss*(Ow - Oss);
    output_GDP    = Yss/GDP;
    NX_nonoil_GDP = NX_nonoil_output     *output_GDP;
    EX_oil_GDP    = EX_oil_output        *output_GDP;
    EX_nonoil_GDP = EX_nonoil_output     *output_GDP;
    IM_nonoil_GDP = IM_nonoil_output     *output_GDP;



    EXoil_EXtotal    = EX_oil_GDP / ( EX_oil_GDP + EX_nonoil_GDP)
    %Trade_GDP        = 2*(EX_oil_GDP + EX_nonoil_GDP)
    OilGDP_totGDP    = ess*POstarss*Ow / GDP
    %ess*PFstarss*CFss/(ess*PFstarss*CFss + CHss )


    %% ---------- FIXED-POINT ON volat ----------
    volat     = 2.6914e-04;%4.0544e-04;             % initial guess
    volat_dcp = 3.2878e-04 ; 
    volat_fxi = 2.2878e-05;
    dX        = 1.0;
    dX_dcp    = 1.0;
    dX_fxi    = 1.0;

    volattemp = zeros(1,100);
    i = 0;

    while abs(dX) > 1e-5 || abs(dX_dcp) > 1e-5  %|| abs(dX_fxi) > 1e-4 
        paramsnew = struct('Cstar', Cstar, 'Ow', Ow, 'gamma', gamma, 'alpha', alpha, ...
            'epsx', epsx, 'beta', beta, 'omega', omega, 'volat', volat,'volat_dcp',volat_dcp, ...
            'epsh', epsh, 'theta', theta, 'Rss', Rss, 'rhor', rhor, ...
            'eta', eta, 'rhon', rhon, 'rhoa', rhoa, 'rhopos', rhopos, ...
            'rhors', rhors, 'rhopf', rhopf, 'Rstarb', Rstarb, ...
            'gammaPi', gammaPi, 'gammaYY', gammaYY, 'Piss', Piss, ...
            'LOss', LOss, 'CFss', CFss, 'Yss', Yss, ...
            'PHstarss',PHstarss,'CHstarss',CHstarss,'ess',ess,'CHss',CHss, ...
            'Lss',Lss,'rerss',rerss,'Oss',Oss,'sigman',sigman);




        % Save params for Dynare .mod to read
        save('steadystate_OEXRPolicy.mat', 'paramsnew');

        % IMPORTANT: Keep workspace across Dynare calls
        % (make sure OEXRPolicy.mod uses the saved .mat)
        dynare OEXRPolicy noclearall

        % Update fixed-point
        dX    = 1e4*(volat     - oo_.var(1,1));
        dX_dcp= 1e4*(volat_dcp - oo_.var(2,2));
        
        volat     = 0.3*volat     + 0.7*oo_.var(1,1);
        volat_dcp = 0.3*volat_dcp + 0.7*oo_.var(2,2);
        i = i+1;
        volattemp(i) = volat;

       % dynare OEXRPolicy_suboptimalfxi noclearall
        dX_fxi    = 1e4*(volat_fxi     - oo_.var(1,1));
        volat_fxi     = 0.3*volat_fxi     + 0.7*oo_.var(1,1);
        
    end


 dynare OEXRPolicy noclearall  
    
    %% ---------- WELFARE / LOSS CALCULATIONS (store per gamma) ----------
    % Indices (use strmatch as in your code, or switch to strcmp depending on Dynare version)
    CH_pos     = strmatch('CH', M_.endo_names, 'exact');
    CF_pos     = strmatch('CF', M_.endo_names, 'exact'); 
    L_pos      = strmatch('L',  M_.endo_names, 'exact'); 

    CH_til_pos = strmatch('CH_til', M_.endo_names, 'exact');
    L_til_pos  = strmatch('L_til',  M_.endo_names, 'exact');

    CH_peg_pos = strmatch('CH_peg', M_.endo_names, 'exact'); 
    L_peg_pos  = strmatch('L_peg',  M_.endo_names, 'exact'); 

    v_peg_pos  = strmatch('v_peg',  M_.endo_names, 'exact');
    Pi_peg_pos = strmatch('Pi_peg', M_.endo_names, 'exact');

    v_noFXI_pos = strmatch('v_noFXI', M_.endo_names, 'exact');
    Pi_pos      = strmatch('Pi',      M_.endo_names, 'exact');
    u_noFXI_pos = strmatch('u_noFXI', M_.endo_names, 'exact');

    v_pegsub_pos      = strmatch('v_pegsub',      M_.endo_names, 'exact');
    Pi_pegsub_pos     = strmatch('Pi_pegsub',     M_.endo_names, 'exact');
    LO_pegsub_pos     = strmatch('LO_pegsub',     M_.endo_names, 'exact');
    LO_til_pos        = strmatch('LO_til',        M_.endo_names, 'exact');
    PHstar_pegsub_pos = strmatch('PHstar_pegsub', M_.endo_names, 'exact');
    PHstar_til_pos    = strmatch('PHstar_til',    M_.endo_names, 'exact');
    
    v_noFXI_dcp_pos   = strmatch('v_noFXI_dcp',    M_.endo_names, 'exact');
    u_noFXI_dcp_pos   = strmatch('u_noFXI_dcp',    M_.endo_names, 'exact');
    Pi_dcp_pos        = strmatch('Pi_dcp',         M_.endo_names, 'exact');
    LO_dcp_pos        = strmatch('LO_dcp',         M_.endo_names, 'exact');
    PHstar_dcp_pos    = strmatch('PHstar_dcp',     M_.endo_names, 'exact');

    Ld_e_pos          = strmatch('Ld_e',          M_.endo_names, 'exact');
    Delta_lnq_pos     = strmatch('Delta_lnq',     M_.endo_names, 'exact');
    Delta_lnq2_pos    = strmatch('Delta_lnq2',     M_.endo_names, 'exact');
    Delta_lnrer_pos   = strmatch('Delta_lnrer',     M_.endo_names, 'exact');
    Delta_lnrer_dcp_pos= strmatch('Delta_lnrer_dcp',     M_.endo_names, 'exact');
    
    ehat_pos          = strmatch('ehat',          M_.endo_names, 'exact');
    qhat_pos          = strmatch('qhat',          M_.endo_names, 'exact');
   
    Ld_e_dcp_pos      = strmatch('Ld_e_dcp',      M_.endo_names, 'exact');

    % Burn-in and length checks
    burn_in = 1000;
    Tsim = size(oo_.endo_simul, 2);
    if Tsim < burn_in+1000
        warning('Simulation length (%d) may be too short for burn-in=%d.', Tsim, burn_in);
    end


    Ld_e      = 0.01*std(oo_.endo_simul(Ld_e_pos,  burn_in:end));
    Ld_e_dcp  = 0.01*std(oo_.endo_simul(Ld_e_dcp_pos,  burn_in:end));
    Delta_lnq = std(oo_.endo_simul(Delta_lnq_pos,  burn_in:end));
    Delta_lnq2 = std(oo_.endo_simul(Delta_lnq2_pos,  burn_in:end));
    Delta_lnrer = std(oo_.endo_simul(Delta_lnrer_pos,  burn_in:end));
    Delta_lnrer_dcp = std(oo_.endo_simul(Delta_lnrer_dcp_pos,  burn_in:end));
   
    ehat  = std(oo_.endo_simul(ehat_pos,  burn_in:end));
    qhat  = std(oo_.endo_simul(qhat_pos,  burn_in:end));
   
    
    % Loss series
    loss_series_peg = (1-gamma) * ( oo_.endo_simul(v_peg_pos,  burn_in:end).^2 + ...
        theta/((1-theta)*(1-beta*theta)) * epsh*(epsh+1) * oo_.endo_simul(Pi_peg_pos, burn_in:end).^2 );

    loss_series_noFXI = (1-gamma) * ( oo_.endo_simul(v_noFXI_pos, burn_in:end).^2 + ...
        theta/((1-theta)*(1-beta*theta)) * epsh*(epsh+1) * log(oo_.endo_simul(Pi_pos, burn_in:end)).^2 ) + ...
        ( gamma + (1-gamma)*(1-alpha)*alpha*Yss/CHss + (1-gamma)*alpha^2*epsx*CHstarss/CHss ) .* ...
        oo_.endo_simul(u_noFXI_pos, burn_in:end).^2;

    loss_series_noFXI_dcp = (1-gamma) * ( oo_.endo_simul(v_noFXI_dcp_pos, burn_in:end).^2 + ...
        theta/((1-theta)*(1-beta*theta)) * epsh*(epsh+1) * log(oo_.endo_simul(Pi_dcp_pos, burn_in:end)).^2 ) + ...
        gamma   *  oo_.endo_simul(u_noFXI_dcp_pos, burn_in:end).^2 + ...
        alpha*(1-alpha)*(1-gamma) * Yss/CHss .* ( log(oo_.endo_simul(LO_dcp_pos, burn_in:end)) - log(oo_.endo_simul(LO_til_pos, burn_in:end)) ).^2 + ...
        ( (1-gamma)*epsx*CHstarss/CHss )     .* ( log(PHstarss) - log(oo_.endo_simul(PHstar_til_pos, burn_in:end)) ).^2;


    loss_series_pegsub = (1-gamma) * ( oo_.endo_simul(v_pegsub_pos, burn_in:end).^2 + ...
        theta/((1-theta)*(1-beta*theta)) * epsh*(epsh+1) * (oo_.endo_simul(Pi_pegsub_pos, burn_in:end)).^2 ) + ...
        alpha*(1-alpha)*(1-gamma) * Yss/CHss .* ( log(oo_.endo_simul(LO_pegsub_pos, burn_in:end)) - log(oo_.endo_simul(LO_til_pos, burn_in:end)) ).^2 + ...
        ( (1-gamma)*epsx*CHstarss/CHss )     .* ( log(oo_.endo_simul(PHstar_pegsub_pos, burn_in:end)) - log(oo_.endo_simul(PHstar_til_pos, burn_in:end)) ).^2;

    
   % dynare OEXRPolicy_suboptimalfxi noclearall
    loss_series_suboptimalFXI = (1-gamma) * ( oo_.endo_simul(v_noFXI_pos, burn_in:end).^2 + ...
            theta/((1-theta)*(1-beta*theta)) * epsh*(epsh+1) * log(oo_.endo_simul(Pi_pos, burn_in:end)).^2 ) + ...
            ( gamma + (1-gamma)*(1-alpha)*alpha*Yss/CHss + (1-gamma)*alpha^2*epsx*CHstarss/CHss ) .* ...
            oo_.endo_simul(u_noFXI_pos, burn_in:end).^2;


    loss_peg            = 1e4*mean(loss_series_peg);
    loss_noFXI          = 1e4*mean(loss_series_noFXI);
    loss_suboptimalFXI  = 1e4*mean(loss_series_suboptimalFXI);
    loss_pegsub         = 1e4*mean(loss_series_pegsub);
    loss_noFXI_dcp      = 1e4*mean(loss_series_noFXI_dcp); 

    lambda_noFXI         = 1 - exp(-loss_noFXI);
    lambda_suboptimalFXI = 1 - exp(-loss_suboptimalFXI);
    lambda_peg           = 1 - exp(-loss_peg);
    lambda_pegsub        = 1 - exp(-loss_pegsub);
    lambda_noFXI_dcp     = 1 - exp(-loss_noFXI_dcp);
    %% ---------- STORE EVERYTHING ----------
    results(jj).gamma        = gamma;
    results(jj).ss           = ss;
    results(jj).volat        = volat;

    results(jj).loss_peg     = loss_peg;
    results(jj).loss_noFXI   = loss_noFXI;
    results(jj).loss_suboptimalFXI   = loss_suboptimalFXI;
    results(jj).loss_pegsub  = loss_pegsub;
    results(jj).loss_noFXI_dcp  = loss_noFXI_dcp;

    results(jj).lambda_peg   = lambda_peg;
    results(jj).lambda_noFXI = lambda_noFXI;
    results(jj).lambda_suboptimalFXI = lambda_suboptimalFXI;
    results(jj).lambda_pegsub= lambda_pegsub;
    results(jj).lambda_noFXI_dcp = lambda_noFXI_dcp;

    results(jj).CHss         = CHss;
    results(jj).Yss          = Yss;
    results(jj).CHstarss     = CHstarss;
    results(jj).LOss         = LOss;

    results(jj).Delta_lnq       = Delta_lnq;
    results(jj).Delta_lnrer     = Delta_lnrer;
    results(jj).Delta_lnrer_dcp = Delta_lnrer_dcp;
    results(jj).Ld_e            = Ld_e;
    results(jj).Ld_e_dcp        = Ld_e_dcp;
    results(jj).weight = gamma + (1-gamma)*(1-alpha)*alpha*Yss/CHss + (1-gamma)*alpha^2*epsx*CHstarss/CHss;
    results(jj).uwedge = mean(oo_.endo_simul(u_noFXI_pos, burn_in:end).^2);
    results(jj).vwedge = mean(oo_.endo_simul(v_noFXI_pos, burn_in:end).^2);
   
    % (Optional) Clean Dynare artifacts between runs to avoid path clashes
    % dynare_cleanup; % or: rmdir('OEXRPolicy/','s'); delete OEXRPolicy*.m etc. if needed
end

%% ---------- SUMMARY TABLE ACROSS gamma ----------
summary = table( ...
    [results.gamma]', ...
    [results.volat]',[results.lambda_peg]', [results.lambda_noFXI]', ...
    [results.lambda_pegsub]',[results.lambda_noFXI_dcp]', [results.lambda_suboptimalFXI]', ...
    [results.loss_peg]',  [results.loss_noFXI]', ...
    [results.loss_pegsub]',[results.loss_noFXI_dcp]',[results.loss_suboptimalFXI]',  [results.Ld_e]', ...
    [results.Ld_e_dcp]',[results.Delta_lnq]',[results.Delta_lnrer]',[results.Delta_lnrer_dcp]', ...
    'VariableNames', {'gamma','volat','lambda_peg','lambda_noFXI', ...
                      'lambda_pegsub','lambda_noFXI_dcp','lambda_suboptimalFXI', ...
                      'loss_peg','loss_noFXI','loss_pegsub', ...
                      'loss_noFXI_dcp','loss_suboptimalFXI','Ld_e','Ld_e_dcp','Delta_lnq','Delta_lnrer','Delta_lnrer_dcp'});

writetable(summary, 'summary_gamma_welfare_all noise.csv');


if true

figure('Name','Convergance');
plot(volattemp);
%%  ****************   Figure 1-a  **************************************
%%  ****************   Figure 1-a  **************************************
%%  ****************   Figure 1-a  **************************************

figure('Name', 'Welfare Cost Comparison');

orange_color = [1, 0.5, 0.1]; 
plot(gamma_v, summary.lambda_peg*100, 'k-o', 'LineWidth', 2.5, 'DisplayName', 'Peg');
hold on;
plot(gamma_v, summary.lambda_pegsub*100, '--d', 'Color', orange_color, 'LineWidth', 2.5, 'DisplayName', 'Peg + Subsidy');
plot(gamma_v, summary.lambda_noFXI*100, 'b--s', 'LineWidth', 2.5, 'DisplayName', 'No FXI');
plot(gamma_v, summary.lambda_noFXI_dcp*100, 'r-.*', 'LineWidth', 2.5, 'MarkerFaceColor', 'r', 'DisplayName', 'No FXI + DCP');
%plot(summary.gamma, summary.lambda_suboptimalFXI*100, 'g-^', 'LineWidth', 2.5, 'MarkerFaceColor', 'g', 'DisplayName', 'Suboptimal FXI');

plot(gamma_v, zeros(size(summary.gamma)), 'k:', 'LineWidth', 2.5, 'DisplayName', 'Optimal Policy');
% --- Add Labels and Formatting ---
xlabel('Share of Oil in Production (1-\alpha)', 'FontSize', 16);
ylabel('Welfare Cost (\lambda, in % CEV)', 'FontSize', 16);
legend('show', 'Location', 'best', 'FontSize', 12); % This command creates the legend box
grid on;
axis tight;

% Optional: Adjust y-axis to start from a clean point if needed
yl = ylim;
ylim([0 yl(2)]); % Force y-axis to start at 0

set(gcf, 'PaperPositionMode', 'auto');
fig_pos = get(gcf, 'PaperPosition');
set(gcf, 'PaperSize', [fig_pos(3) fig_pos(4)]);
ax = gca;
set(gcf,'Units','inches','Position',[0 0 7 5]);  % width x height
set(gcf,'Color','none');
set(gca, 'FontSize', 14);
set(gca, 'LineWidth', 1.2);
exportgraphics(gcf,'Welfare_gamma.pdf', 'ContentType','vector','BackgroundColor','none');

end






%%  ****************   Figure 1-c  **************************************
%%  ****************   Figure 1-c  **************************************
%%  ****************   Figure 1-c  **************************************

    data_noise = readtable('summary_gamma_welfare_all noise.csv');
    data_no_noise = readtable('summary_gamma_welfare_all NOnoise.csv');

    lambda_noFXI_noise = data_noise.lambda_noFXI;
    lambda_noFXI_dcp_noise = data_noise.lambda_noFXI_dcp;

    lambda_noFXI_no_noise = data_no_noise.lambda_noFXI;
    lambda_noFXI_dcp_no_noise = data_no_noise.lambda_noFXI_dcp;

my_lambda_noFXI = (lambda_noFXI_noise - lambda_noFXI_no_noise) * 100;
my_lambda_noFXI_dcp = (lambda_noFXI_dcp_noise - lambda_noFXI_dcp_no_noise) * 100;

% For the peg regimes, the difference should be zero, as noise has no effect.
my_lambda_peg = zeros(size(gamma_v));
my_lambda_pegsub = zeros(size(gamma_v));

figure('Name', 'Welfare Cost Comparison');

plot(gamma_v, my_lambda_peg, 'k--d', 'LineWidth', 2.5, 'MarkerFaceColor', [0.47 0.67 0.19], 'DisplayName', 'Peg');
hold on;
plot(gamma_v, my_lambda_pegsub, '-.*', 'Color', [0.9, 0.5, 0.1], 'LineWidth', 2.5, 'DisplayName', 'Peg + Subsidy'); % Changed style for clarity
plot(gamma_v, my_lambda_noFXI, 'b-o', 'LineWidth', 2.5, 'DisplayName', 'No FXI'); % Changed style for clarity
plot(gamma_v, my_lambda_noFXI_dcp, 'r--s', 'LineWidth', 2.5, 'DisplayName', 'No FXI + DCP');
% --- Add Labels and Formatting ---
xlabel('Share of Oil in Production (1-\alpha)', 'FontSize', 16);
ylabel('Welfare Cost (% CEV)', 'FontSize', 16);
legend('show', 'Location', 'best', 'FontSize', 12); % This command creates the legend box
grid on;
axis tight;

% Optional: Adjust y-axis to start from a clean point if needed
yl = ylim;
ylim([0 yl(2)]); % Force y-axis to start at 0

set(gcf, 'PaperPositionMode', 'auto');
fig_pos = get(gcf, 'PaperPosition');
set(gcf, 'PaperSize', [fig_pos(3) fig_pos(4)]);
ax = gca;
set(gcf,'Units','inches','Position',[0 0 7 5]);  % width x height
set(gcf,'Color','none');
set(gca, 'FontSize', 14);
set(gca, 'LineWidth', 1.2);
exportgraphics(gcf,'Welfare_gamma_NOISECOMPARISON.pdf','ContentType','vector','BackgroundColor','none');